Enrichment x stress MA
Supplementary Material
Setting-up
Loading packages
#devtools::install_github('Mikata-Project/ggthemr', force = TRUE)
pacman::p_load(tidyverse,
here,
metafor,
clubSandwich,
orchaRd,
MuMIn,
patchwork,
GoodmanKruskal,
networkD3,
ggplot2,
visdat,
ggalluvial,
ggthemr,
cowplot,
grDevices,
png,
grid,
gridGraphics,
pander)##
## The downloaded binary packages are in
## /var/folders/b0/5_1bfml50hj4c3_9l0kzfzyr0000gn/T//RtmpdnAq06/downloaded_packages
# needed for model selection using MuMIn within metafor
eval(metafor:::.MuMIn)Loading data and functions
This loads the unprocessed datafile and custom functions including: calculating ‘focal’ and ‘pair-wise’ effect sizes and variance calculating SMD *creating VCV of variance
dat <- read_csv(here("Data","Data_raw.csv"))
# Load custom function to extract data
source(here("R/Functions.R")) Data organisation
removing study (Wang et al_2020) with negative values as lnRR cannot be calculated with negative values rounding down sample sizes that are reported as decimals due to averaging n across treatments (rounding down is more conservative), getting effect sizes from function ’flipping’ effect sizes so that all effect sizes are higher values = individuals do better and learning/memory *assigning human readable terms, and creating VCV of variance
# removing study with negative values as these are unable to be used for lnRR
dat <- droplevels(dat[!dat$First_author == 'Wang',])
#rounding down sample sizes
dat$CC_n <- floor(dat$CC_n)
dat$EC_n <- floor(dat$EC_n)
dat$CS_n <- floor(dat$CS_n)
dat$ES_n <- floor(dat$CS_n)
# 'Focal' effect_size
effect_size <- with(dat, mapply(effect_set,
CC_n ,
CC_mean,
CC_SD,
EC_n,
EC_mean,
EC_SD,
CS_n,
CS_mean,
CS_SD,
ES_n,
ES_mean,
ES_SD,
percent = Response_percent,
SIMPLIFY = FALSE))
effect_size <- map_dfr(effect_size, I)
# 'Pairwise' effect size
effect_size2 <- with(dat, mapply(effect_set2,
CC_n ,
CC_mean,
CC_SD,
EC_n,
EC_mean,
EC_SD,
CS_n,
CS_mean,
CS_SD,
ES_n,
ES_mean,
ES_SD,
percent = Response_percent,
SIMPLIFY = FALSE))
effect_size2 <- map_dfr(effect_size2, I)
full_info <- which(complete.cases(effect_size) == TRUE)
# adding effect sizes as column
dat <- bind_cols(dat, effect_size, effect_size2)
dat <- dat[full_info, ]
#Flipping 'lower is better' to 'higher is better' effect sizes
#flipping lnRR for values where higher = worse
dat$lnRR_Ea <- ifelse(dat$Response_direction == 2, dat$lnRR_E*-1,ifelse(is.na(dat$Response_direction) == TRUE, NA, dat$lnRR_E))
# currently NAswhich causes error
dat$lnRR_Sa <- ifelse(dat$Response_direction == 2, dat$lnRR_S*-1,ifelse(is.na(dat$Response_direction) == TRUE, NA, dat$lnRR_S)) # currently NAswhich causes error
dat$lnRR_ESa <- ifelse(dat$Response_direction == 2, dat$lnRR_ES*-1,ifelse(is.na(dat$Response_direction) == TRUE, NA, dat$lnRR_ES)) # currently NAswhich causes error
#flipping 'pure effect sizes'
dat$lnRR_E2a <- ifelse(dat$Response_direction == 2, dat$lnRR_E2*-1,ifelse(is.na(dat$Response_direction) == TRUE, NA, dat$lnRR_E2)) # currently NAswhich causes error
dat$lnRR_S2a <- ifelse(dat$Response_direction == 2, dat$lnRR_S2*-1,ifelse(is.na(dat$Response_direction) == TRUE, NA, dat$lnRR_S2)) # currently NAswhich causes error
dat$lnRR_ES2a <- ifelse(dat$Response_direction == 2, dat$lnRR_ES2*-1,ifelse(is.na(dat$Response_direction) == TRUE, NA, dat$lnRR_ES2)) # currently NAswhich causes error
dat$lnRR_E3a <- ifelse(dat$Response_direction == 2, dat$lnRR_E3*-1,ifelse(is.na(dat$Response_direction) == TRUE, NA, dat$lnRR_E3)) # currently NAswhich causes error
dat$lnRR_S3a <- ifelse(dat$Response_direction == 2, dat$lnRR_S3*-1,ifelse(is.na(dat$Response_direction) == TRUE, NA, dat$lnRR_S3)) # currently NAswhich causes error
#flipping SMD
dat$SMD_Ea <- ifelse(dat$Response_direction == 2, dat$SMD_E*-1,ifelse(is.na(dat$Response_direction) == TRUE, NA, dat$SMD_E)) # currently NAswhich causes error
dat$SMD_Sa <- ifelse(dat$Response_direction == 2, dat$SMD_S*-1,ifelse(is.na(dat$Response_direction) == TRUE, NA, dat$SMD_S)) # currently NAswhich causes error
dat$SMD_ESa <- ifelse(dat$Response_direction == 2, dat$SMD_ES*-1,ifelse(is.na(dat$Response_direction) == TRUE, NA, dat$SMD_ES))
# assigning human readable terms
dat <- dat %>% mutate(Type_assay = case_when(Type_assay == 1 ~ "Habituation",
Type_assay == 2 ~ "Conditioning",
Type_assay == 3 ~ "Recognition",
Type_assay == 4 ~ "Unclear"),
Learning_vs_memory = case_when(Learning_vs_memory == 1 ~ "Learning",
Learning_vs_memory == 2 ~ "Memory",
Learning_vs_memory == 3 ~ "Habituation"),
Type_reinforcement = case_when(Type_reinforcement== 1 ~"Appetitive",
Type_reinforcement== 2 ~ "Aversive",
Type_reinforcement== 3 ~ "Not applicable",
Type_reinforcement== 4 ~ "Unclear"),
Type_stress_exposure = case_when(Type_stress_exposure == 1 ~ "Density",
Type_stress_exposure == 2 ~ "Scent",
Type_stress_exposure == 3 ~ "Shock",
Type_stress_exposure == 4 ~ "Exertion",
Type_stress_exposure == 5 ~ "Restraint",
Type_stress_exposure == 6 ~ "MS",
Type_stress_exposure == 7 ~ "Circadian rhythm",
Type_stress_exposure == 8 ~ "Noise",
Type_stress_exposure == 9 ~ "Other",
Type_stress_exposure == 10 ~ "Combination",
Type_stress_exposure == 11 ~ "unclear"),
Age_stress_exposure = case_when(Age_stress_exposure == 1 ~ "Prenatal",
Age_stress_exposure == 2 ~ "Early postnatal",
Age_stress_exposure == 3 ~ "Adolescent",
Age_stress_exposure == 4 ~ "Adult",
Age_stress_exposure == 5 ~ "Unclear"),
Stress_duration = case_when(Stress_duration == 1 ~ "Acute",
Stress_duration == 2 ~ "Chronic",
Stress_duration == 3 ~ "Intermittent",
Stress_duration == 4 ~ "Unclear"),
EE_social = case_when(EE_social == 1 ~ "Social",
EE_social== 2 ~ "Non-social",
EE_social == 3 ~ "Unclear"),
EE_exercise = case_when(EE_exercise == 1 ~ "Exercise",
EE_exercise == 2 ~ "No exercise"),
Age_EE_exposure = case_when(Age_EE_exposure == 1 ~ "Prenatal",
Age_EE_exposure == 2 ~ "Early postnatal",
Age_EE_exposure == 3 ~ "Adolescent",
Age_EE_exposure == 4 ~ "Adult",
Age_EE_exposure == 5 ~ "Unclear"),
Exposure_order = case_when(Exposure_order == 1 ~ "Stress first",
Exposure_order == 2 ~ "Enrichment first",
Exposure_order == 3 ~ "Concurrently",
Exposure_order == 4 ~ "Unclear"),
Age_assay = case_when(Age_assay == 1 ~ "Early postnatal",
Age_assay == 2 ~ "Adolescent",
Age_assay == 3 ~ "Adult",
Age_assay == 4 ~ "Unclear"),
Sex = case_when(Sex == 1 ~ "Female",
Sex == 2 ~ "Male",
Sex == 3 ~ "Mixed",
Sex == 4 ~ "Unclear"),
Type_EE_exposure = case_when(Type_EE_exposure == 1 ~ "Nesting material",
Type_EE_exposure == 2 ~ "Objects",
Type_EE_exposure == 3 ~ "Cage complexity",
Type_EE_exposure == 4 ~ "Wheel/trademill",
Type_EE_exposure == 5 ~ "Combination",
Type_EE_exposure == 6 ~ "Other",
Type_EE_exposure == 7 ~ "Unclear"),
ROB_blinding = case_when(ROB_blinding == 1 ~ "Yes",
ROB_blinding == 2 ~ "No",
ROB_blinding == 3 ~ "Unclear"),
ROB_randomisation = case_when(ROB_randomisation == 1 ~ "Yes",
ROB_randomisation == 2 ~ "No",
ROB_randomisation == 3 ~ "Unclear"))
#making variance VCVs
VCV_E <- impute_covariance_matrix(vi = dat$lnRRV_E, cluster = dat$Study_ID, r = 0.5)
VCV_S <- impute_covariance_matrix(vi = dat$lnRRV_S, cluster = dat$Study_ID, r = 0.5)
VCV_ES <- impute_covariance_matrix(vi = dat$lnRRV_ES, cluster = dat$Study_ID, r = 0.5)
VCV_Ea <- impute_covariance_matrix(vi = dat$SMDV_E, cluster = dat$Study_ID, r = 0.5)
VCV_Sa <- impute_covariance_matrix(vi = dat$SMDV_S, cluster = dat$Study_ID, r = 0.5)
VCV_ESa <- impute_covariance_matrix(vi = dat$SMDV_ES, cluster = dat$Study_ID, r = 0.5)
#write.csv(dat, file = here("Data", 'Data_processed.csv'), row.names = TRUE)Data exploration
General
number of effect sizes number of studies *publication year range
#Number of effect sizes
length(unique(dat$ES_ID))
#Number of studies
length(unique(dat$Study_ID))
#Publication years
min(dat$Year_published)
max(dat$Year_published)Visual of missing data
plot_missing <- vis_miss(dat) +
theme(plot.title = element_text(hjust = 0.5, vjust = 3),
plot.margin = margin(t = 0.5, r = 3, b = 1, l = 1, unit = "cm")) +
ggtitle("Missing data in the selected predictors") #no missing values
plot_missing#useGoodman and Kruskal’s τ measure of association between categorical predictor variables (function from package GoodmanKruskal: https://cran.r-project.org/web/packages/GoodmanKruskal/vignettes/GoodmanKruskal.html)
#GKmatrix <- GKtauDataframe(subset(dat, select = c("Sex", "Type_assay", "Learning_vs_memory", #"Type_reinforcement", "Type_stress_exposure", "Age_stress_exposure", "Stress_duration", #"EE_social_HR", "EE_exercise", "Age_EE_exposure", "Exposure_order", "Age_assay")))
#plot(GKmatrix)
#simple pairwise contingency tables
# table(dat$Type_assay, dat$Type_reinforcement)
# table(dat$Age_stress_exposure, dat$Age_EE_exposure)
# table(dat$Type_stress_exposure, dat$Age_stress_exposure)
# table(dat$Type_stress_exposure, dat$Age_assay)
# table(dat$Type_stress_exposure, dat$Stress_duration)Alluvial diagrams
Shows the relationship/nestedness of different data elements (used for produce Fig. 2)
Subjects info: species-strain-sex
freq_A <- as.data.frame(table(dat$Sex, dat$Common_species, dat$Strain)) %>% rename(Sex = Var1, Species = Var2, Strain = Var3) #make a data frame of frequencies for three selected variables
is_alluvia_form(as.data.frame(freq_A), axes = 1:3, silent = TRUE)
p1 <- ggplot(data = freq_A,
aes(axis1 = Sex, axis2 = Species, axis3 = Strain, y = Freq)) +
geom_alluvium(aes(fill = Sex)) +
geom_flow() +
geom_stratum(aes(fill = Sex)) +
geom_text(stat = "stratum", aes(label = after_stat(stratum))) +
#theme_minimal() +
theme_void() +
theme(legend.position = "none",
plot.title = element_text(hjust = 0, vjust = 3),
axis.title.x = element_text(),
axis.text.x = element_text(face="bold"),
plot.margin = unit(c(1, 1, 0, 1), "cm")) +
scale_x_discrete(limits = c("Sex", "Species", "Strain"), expand = c(0.15, 0.05), position = "top") +
ggtitle("A study subjects")
p1Stress info: age-duration-type stress
freq_C <- as.data.frame(table(dat$Age_stress_exposure, dat$Stress_duration, dat$Type_stress_exposure)) %>% rename(Age_stress = Var1, Duration_stress = Var2, Type_stress = Var3) #make a data frame of frequencies for three selected variables
is_alluvia_form(as.data.frame(freq_C), axes = 1:3, silent = TRUE)
p3 <- ggplot(data = freq_C,
aes(axis1 = Age_stress, axis2 = Duration_stress, axis3 = Type_stress, y = Freq)) +
geom_alluvium(aes(fill = Age_stress)) +
geom_flow() +
geom_stratum(aes(fill = Age_stress)) +
geom_text(stat = "stratum", aes(label = after_stat(stratum))) +
#theme_minimal() +
theme_void() +
theme(legend.position = "none",
plot.title = element_text(hjust = 0, vjust = 3),
axis.title.x = element_text(),
axis.text.x = element_text(face="bold"),
plot.margin = unit(c(1, 1, 0, 1), "cm")) +
scale_x_discrete(limits = c("Age", "Duration", "Type"), expand = c(0.1, 0.1), position = "top") +
ggtitle("C stress exposure")
p3Assay info: LvsM-type-reinforcement
freq_D <- as.data.frame(table(dat$Learning_vs_memory, dat$Type_assay, dat$Type_reinforcement)) %>% rename(Learning_Memory = Var1, Type = Var2, Reinforcement = Var3) #make a data frame of frequencies for three selected variables
is_alluvia_form(as.data.frame(freq_D), axes = 1:3, silent = TRUE)
p4 <- ggplot(data = freq_D,
aes(axis1 = Learning_Memory, axis2 = Type, axis3 = Reinforcement, y = Freq)) +
geom_alluvium(aes(fill = Learning_Memory)) +
geom_flow() +
geom_stratum(aes(fill = Learning_Memory)) +
geom_text(stat = "stratum", aes(label = after_stat(stratum))) +
#theme_minimal() +
theme_void() +
theme(legend.position = "none",
plot.title = element_text(hjust = 0, vjust = 3),
axis.title.x = element_text(),
axis.text.x = element_text(face="bold"),
plot.margin = unit(c(1, 1, 0, 1), "cm")) +
scale_x_discrete(limits = c("Learning_Memory", "Type", "Reinforcement"), expand = c(0.1, 0.1), position = "top") +
ggtitle("D cognitive assay")
p4Combined plot
Used for Fig. 2
#p1 + scale_fill_brewer(palette = "Set3") #Pastel1
(p1 + scale_fill_brewer(palette = "Set3")) / (p2 + scale_fill_brewer(palette = "Set3")) / (p3 + scale_fill_brewer(palette = "Set3")) / (p4 + scale_fill_brewer(palette = "Set3")) + plot_layout(ncol = 1, heights = c(1,1,1,1,1))#ggsave(file = "./figs/Alluvial_diagrams_v2.pdf", width = 10, height = 12, units = "cm", dpi = 300, scale = 2) #, device = cairo_pdf)Modelling with lnRR
Environmental enrichment
Meta-analysis
#dat <- read_csv(here("Data","Data_processed.csv"))
mod_E0 <- rma.mv(yi = lnRR_Ea, V = VCV_E, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat)
summary(mod_E0) ##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -9.7836 19.5672 27.5672 37.6106 28.0323
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0082 0.0905 30 no Study_ID
## sigma^2.2 0.0345 0.1858 92 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Heterogeneity:
## Q(df = 91) = 809.2712, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## 0.1860 0.0320 5.8116 91 <.0001 0.1224 0.2496 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
i2_ml(mod_E0) ## I2_total I2_Study_ID I2_ES_ID I2_Strain
## 9.310673e-01 1.786605e-01 7.524068e-01 6.097041e-10
orchard_plot(mod_E0, mod = "Int", xlab = "lnRR", alpha=0.4) + # Orchard plot
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5)+ # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2)+ # confidence intervals
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_colour_manual(values = "darkorange")+ # change colours
scale_fill_manual(values="darkorange")+
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 24), # change font sizes
legend.title = element_text(size = 15),
legend.text = element_text(size = 13)) Meta-regression: uni-moderator
Learning vs Memory
Was the response learning or memory?
mod_E1 <- rma.mv(yi = lnRR_Ea, V = VCV_E, mod = ~Learning_vs_memory-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat)
summary(mod_E1) ##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -9.3793 18.7586 28.7586 41.2577 29.4729
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0076 0.0873 30 no Study_ID
## sigma^2.2 0.0340 0.1843 92 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 90) = 802.5794, p-val < .0001
##
## Test of Moderators (coefficients 1:2):
## F(df1 = 2, df2 = 90) = 18.2861, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb
## Learning_vs_memoryLearning 0.2303 0.0450 5.1227 90 <.0001 0.1410
## Learning_vs_memoryMemory 0.1624 0.0355 4.5713 90 <.0001 0.0919
## ci.ub
## Learning_vs_memoryLearning 0.3197 ***
## Learning_vs_memoryMemory 0.2330 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_E1) ## R2_marginal R2_coditional
## 0.02572583 1.00000000
LvsM_E<- orchard_plot(mod_E1, mod = "Learning_vs_memory", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
xlim(-0.5, 2) +
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 15), # change font sizes
legend.title = element_text(size = 10),
legend.text = element_text(size = 10))
LvsM_EType of assay
The broad category of the type of assay used to measure learning or memory
dat1 <- filter(dat, Type_assay %in% c("Recognition", "Habituation", "Conditioning"))
VCV_E1 <- impute_covariance_matrix(vi = dat1$lnRRV_E, cluster = dat$Study_ID, r = 0.5)
mod_E2 <- rma.mv(yi = lnRR_Ea, V = VCV_E1, mod = ~Type_assay-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat1)
summary(mod_E2)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -7.8496 15.6993 27.6993 42.6311 28.7237
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0099 0.0995 30 no Study_ID
## sigma^2.2 0.0321 0.1792 92 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 89) = 661.8903, p-val < .0001
##
## Test of Moderators (coefficients 1:3):
## F(df1 = 3, df2 = 89) = 12.7993, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## Type_assayConditioning 0.2167 0.0351 6.1650 89 <.0001 0.1468 0.2865
## Type_assayHabituation 0.1821 0.1147 1.5886 89 0.1157 -0.0457 0.4100
## Type_assayRecognition 0.0554 0.0640 0.8659 89 0.3889 -0.0717 0.1826
##
## Type_assayConditioning ***
## Type_assayHabituation
## Type_assayRecognition
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_E2) ## R2_marginal R2_coditional
## 0.07376134 1.00000000
Learning_E <- orchard_plot(mod_E2, mod = "Type_assay", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
xlim(-0.5, 2) +
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 15), # change font sizes
legend.title = element_text(size = 10),
legend.text = element_text(size = 10))
Learning_EType of reinforcement
If conditioning was used, was aversive or appetitive reinforcement used?
dat2 <- filter(dat, Type_reinforcement %in% c("Appetitive", "Aversive", "Not applicable"))
VCV_E2 <- impute_covariance_matrix(vi = dat2$lnRRV_E, cluster = dat2$Study_ID, r = 0.5)
mod_E3 <- rma.mv(yi = lnRR_Ea, V = VCV_E2, mod = ~ Type_reinforcement-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat2)
summary(mod_E3)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -7.8702 15.7405 27.7405 42.6723 28.7649
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0109 0.1046 30 no Study_ID
## sigma^2.2 0.0319 0.1787 92 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 89) = 764.5984, p-val < .0001
##
## Test of Moderators (coefficients 1:3):
## F(df1 = 3, df2 = 89) = 12.4138, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb
## Type_reinforcementAppetitive 0.2175 0.0726 2.9942 89 0.0036 0.0732
## Type_reinforcementAversive 0.2191 0.0410 5.3456 89 <.0001 0.1377
## Type_reinforcementNot applicable 0.0812 0.0560 1.4504 89 0.1505 -0.0300
## ci.ub
## Type_reinforcementAppetitive 0.3618 **
## Type_reinforcementAversive 0.3006 ***
## Type_reinforcementNot applicable 0.1924
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_E3) ## R2_marginal R2_coditional
## 0.07720103 1.00000000
Reinforcement_E <-orchard_plot(mod_E3, mod = "Type_reinforcement", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
xlim(-0.5, 2) +
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 15), # change font sizes
legend.title = element_text(size = 10),
legend.text = element_text(size = 10))
Reinforcement_EAge of enrichment
The age when individuals were exposed to environmental enrichment
dat3 <- filter(dat, Age_EE_exposure %in% c("Adult", "Adolescent"))
VCV_E3 <- impute_covariance_matrix(vi = dat3$lnRRV_E, cluster = dat3$Study_ID, r = 0.5)
mod_E4 <- rma.mv(yi = lnRR_Ea, V = VCV_E3, mod = ~Age_EE_exposure-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat3)
summary(mod_E4) ##
## Multivariate Meta-Analysis Model (k = 88; method: REML)
##
## logLik Deviance AIC BIC AICc
## -6.9490 13.8980 23.8980 36.1698 24.6480
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0078 0.0885 29 no Study_ID
## sigma^2.2 0.0327 0.1808 88 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 86) = 782.6092, p-val < .0001
##
## Test of Moderators (coefficients 1:2):
## F(df1 = 2, df2 = 86) = 18.9060, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## Age_EE_exposureAdolescent 0.1799 0.0370 4.8553 86 <.0001 0.1062 0.2535
## Age_EE_exposureAdult 0.2340 0.0620 3.7734 86 0.0003 0.1107 0.3573
##
## Age_EE_exposureAdolescent ***
## Age_EE_exposureAdult ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_E4) ## R2_marginal R2_coditional
## 0.01127347 1.00000000
Age_E<- orchard_plot(mod_E4, mod = "Age_EE_exposure", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
xlim(-0.5, 2) +
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 15), # change font sizes
legend.title = element_text(size = 10),
legend.text = element_text(size = 10))
Age_EExercise enrichment
Did enrichment involve the addition of apparatus for voluntary exercise (e.g., a wheel or treadmill)?
mod_E5<- rma.mv(yi = lnRR_Ea, V = VCV_E, mod = ~EE_exercise-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat)
summary(mod_E5)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -10.0303 20.0606 30.0606 42.5597 30.7749
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0096 0.0980 30 no Study_ID
## sigma^2.2 0.0345 0.1857 92 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 90) = 807.7169, p-val < .0001
##
## Test of Moderators (coefficients 1:2):
## F(df1 = 2, df2 = 90) = 16.1666, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## EE_exerciseExercise 0.1849 0.0407 4.5464 90 <.0001 0.1041 0.2657
## EE_exerciseNo exercise 0.1900 0.0556 3.4151 90 0.0010 0.0795 0.3005
##
## EE_exerciseExercise ***
## EE_exerciseNo exercise ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_E5) ## R2_marginal R2_coditional
## 0.0001360993 0.9999999987
Exercise_E <-orchard_plot(mod_E5, mod = "EE_exercise", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
xlim(-0.5, 2) +
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 15), # change font sizes
legend.title = element_text(size = 10),
legend.text = element_text(size = 10))
Exercise_EMulti-moderator model
The akaike weights for the top set of models with AIC < 6
# filter data so that all K < 5 are removed
dat_Efm <- dat %>%
filter(Type_assay %in% c("Recognition", "Habituation", "Conditioning"),
Type_reinforcement %in% c("Appetitive", "Aversive", "Not applicable"),
EE_social %in% c("Social", "Non-social"),
Age_EE_exposure %in% c("Adult", "Adolescent"))
VCV_Efm <- impute_covariance_matrix(vi = dat_Efm$lnRRV_E, cluster = dat_Efm$Study_ID, r = 0.5)
mod_Efm <- rma.mv(yi = lnRR_Sa, V = VCV_Efm, mod = ~Type_assay-1 + Learning_vs_memory + Type_reinforcement + EE_social + EE_exercise + Age_EE_exposure , random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat_Efm)
#summary(mod_Efm)
#r2_ml(mod_Efm)
res_Efm <- dredge(mod_Efm, trace=2)
saveRDS(res_Efm, file = here("Rdata", "res_Efm.rds"))
# also saving the full model and data
saveRDS(mod_Efm, file = here("Rdata", "mod_Efm.rds"))
saveRDS(dat_Efm, file = here("Rdata", "dat_Efm.rds"))dat_Efm <- readRDS(file = here("Rdata", "dat_Efm.rds"))
mod_Efm <- readRDS(file = here("Rdata", "mod_Efm.rds"))
res_Efm <- readRDS(file = here("Rdata", "res_Efm.rds"))
res_Efm2<- subset(res_Efm, delta <= 6, recalc.weights=FALSE)
importance(res_Efm2)## Age_EE_exposure Type_assay EE_exercise EE_social
## Sum of weights: 0.63 0.36 0.16 0.14
## N containing models: 11 10 5 5
## Learning_vs_memory Type_reinforcement
## Sum of weights: 0.10 0.07
## N containing models: 4 4
Publication bias & sensitivity analysis
Publication bias tests
# funnel plot
Funnel_E<-funnel(mod_Efm, xlab = "lnRR", ylab = "Standard Error")Funnel_E
#year published was scaled previously under stress PB
dat_Efm$sqrt_inv_e_n <- with(dat_Efm, sqrt(1/CC_n + 1/EC_n + 1/ES_n + 1/CS_n))
PB_MR_E<- rma.mv(lnRR_Sa, lnRRV_S, mods = ~1 + sqrt_inv_e_n + Learning_vs_memory + Year_published + Type_assay + Type_reinforcement + EE_social + EE_exercise + Age_stress_exposure, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain), method = "REML", test = "t",
data = dat_Efm)
estimates_PB_MR_E<- estimates.CI(PB_MR_E)
#estimates_PB_MR_ELeave-one-out analysis
dat$Study_ID <- as.factor(dat$Study_ID)
LeaveOneOut_effectsize <- list()
for(i in 1:length(levels(dat$Study_ID))){
LeaveOneOut_effectsize[[i]] <- rma.mv(yi = lnRR_Ea, V = lnRRV_E,
random = list(~1 | Study_ID,~1| ES_ID, ~1 | Strain),
method = "REML", data = dat[dat$Study_ID
!= levels(dat$Study_ID)[i], ])}
# writing function for extracting est, ci.lb, and ci.ub from all models
est.func <- function(mod_E0){
df <- data.frame(est = mod_E0$b, lower = mod_E0$ci.lb, upper = mod_E0$ci.ub)
return(df)
}
#using dplyr to form data frame
MA_CVR_E <- lapply(LeaveOneOut_effectsize, function(x) est.func(x))%>% bind_rows %>% mutate(left_out = levels(dat$Study_ID))
saveRDS(MA_CVR_E,file = here("Rdata", "MA_CVR_E.rds"))#telling ggplot to stop reordering factors
MA_CVR_E <- readRDS(file = here("Rdata", "MA_CVR_E.rds"))
MA_CVR_E$left_out<- as.factor(MA_CVR_E$left_out)
MA_CVR_E$left_out<-factor(MA_CVR_E$left_out, levels = MA_CVR_E$left_out)
#plotting
leaveoneout_E <- ggplot(MA_CVR_E) +
geom_hline(yintercept = 0, lty = 2, lwd = 1) +
geom_hline(yintercept = mod_E0$ci.lb, lty = 3, lwd = 0.75, colour = "black") +
geom_hline(yintercept = mod_E0$b, lty = 1, lwd = 0.75, colour = "black") +
geom_hline(yintercept = mod_E0$ci.ub, lty = 3, lwd = 0.75, colour = "black") +
geom_pointrange(aes(x = left_out, y = est, ymin = lower, ymax = upper)) +
xlab("Study left out") +
ylab("lnRR, 95% CI") +
coord_flip() +
theme(panel.grid.minor = element_blank())+
theme_bw() + theme(panel.grid.major = element_blank()) +
theme(panel.grid.minor.x = element_blank() ) +
theme(axis.text.y = element_text(size = 6))
leaveoneout_Edat$Study_ID <- as.integer(dat$Study_ID)Stress
Meta-analysis
mod_S0 <- rma.mv(yi = lnRR_Sa, V = VCV_S, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat)
summary(mod_S0) ##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -14.1156 28.2312 36.2312 46.2747 36.6964
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0099 0.0993 30 no Study_ID
## sigma^2.2 0.0377 0.1941 92 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Heterogeneity:
## Q(df = 91) = 946.9234, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## -0.1052 0.0337 -3.1217 91 0.0024 -0.1721 -0.0383 **
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
i2_ml(mod_S0) ## I2_total I2_Study_ID I2_ES_ID I2_Strain
## 9.376124e-01 1.946168e-01 7.429955e-01 1.975455e-10
orchard_plot(mod_S0, mod = "Int", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 24), # change font sizes
legend.title = element_text(size = 15),
legend.text = element_text(size = 13)) Meta-regression: uni-moderator
Learning vs Memory
Was the response learning or memory?
mod_S1 <- rma.mv(yi = lnRR_Sa, V = VCV_S, mod = ~Learning_vs_memory-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat)
summary(mod_S1) ##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -14.5281 29.0562 39.0562 51.5553 39.7705
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0094 0.0970 30 no Study_ID
## sigma^2.2 0.0388 0.1969 92 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 90) = 946.3930, p-val < .0001
##
## Test of Moderators (coefficients 1:2):
## F(df1 = 2, df2 = 90) = 5.0145, p-val = 0.0086
##
## Model Results:
##
## estimate se tval df pval ci.lb
## Learning_vs_memoryLearning -0.1211 0.0476 -2.5423 90 0.0127 -0.2157
## Learning_vs_memoryMemory -0.0974 0.0380 -2.5648 90 0.0120 -0.1728
## ci.ub
## Learning_vs_memoryLearning -0.0265 *
## Learning_vs_memoryMemory -0.0219 *
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_S1) ## R2_marginal R2_coditional
## 0.002776034 1.000000000
LvsM_S <- orchard_plot(mod_S1, mod = "Learning_vs_memory", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 15), # change font sizes
legend.title = element_text(size = 10),
legend.text = element_text(size = 10))
LvsM_S Type of assay
The broad category of the type of assay used to measure learning or memory
dat$Type_assay<-as.factor(dat$Type_assay)
VCV_S1 <- impute_covariance_matrix(vi = dat1$lnRRV_S, cluster = dat$Study_ID, r = 0.5)
mod_S2 <- rma.mv(yi = lnRR_Sa, V = VCV_S1, mod = ~Type_assay-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat1)
summary(mod_S2)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -9.8028 19.6055 31.6055 46.5374 32.6299
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0161 0.1271 30 no Study_ID
## sigma^2.2 0.0279 0.1671 92 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 89) = 723.4973, p-val < .0001
##
## Test of Moderators (coefficients 1:3):
## F(df1 = 3, df2 = 89) = 6.7053, p-val = 0.0004
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## Type_assayConditioning -0.0981 0.0375 -2.6192 89 0.0104 -0.1725 -0.0237
## Type_assayHabituation -0.4615 0.1126 -4.0969 89 <.0001 -0.6853 -0.2377
## Type_assayRecognition -0.0534 0.0645 -0.8287 89 0.4095 -0.1816 0.0747
##
## Type_assayConditioning *
## Type_assayHabituation ***
## Type_assayRecognition
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_S2) ## R2_marginal R2_coditional
## 0.1853359 1.0000000
Learning_S <-orchard_plot(mod_S2, mod = "Type_assay", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 15), # change font sizes
legend.title = element_text(size = 10),
legend.text = element_text(size = 10))
Learning_SType of reinforcement
If conditioning was used, was aversive or appetitive reinforcement used?
VCV_S2 <- impute_covariance_matrix(vi = dat2$lnRRV_S, cluster = dat$Study_ID, r = 0.5)
mod_S3 <- rma.mv(yi = lnRR_Sa, V = VCV_S2, mod = ~ Type_reinforcement-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat2)
summary(mod_S3)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -13.8810 27.7621 39.7621 54.6939 40.7864
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0103 0.1016 30 no Study_ID
## sigma^2.2 0.0387 0.1966 92 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 89) = 920.8439, p-val < .0001
##
## Test of Moderators (coefficients 1:3):
## F(df1 = 3, df2 = 89) = 3.7942, p-val = 0.0130
##
## Model Results:
##
## estimate se tval df pval
## Type_reinforcementAppetitive -0.1846 0.0749 -2.4649 89 0.0156
## Type_reinforcementAversive -0.0730 0.0427 -1.7081 89 0.0911
## Type_reinforcementNot applicable -0.1172 0.0590 -1.9851 89 0.0502
## ci.lb ci.ub
## Type_reinforcementAppetitive -0.3334 -0.0358 *
## Type_reinforcementAversive -0.1579 0.0119 .
## Type_reinforcementNot applicable -0.2345 0.0001 .
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_S3) ## R2_marginal R2_coditional
## 0.0366428 1.0000000
Reinforcement_S <-orchard_plot(mod_S3, mod = "Type_reinforcement", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 15), # change font sizes
legend.title = element_text(size = 10),
legend.text = element_text(size = 10))
Reinforcement_SAge of stress
The age when individuals were exposed to stress
mod_S4 <-rma.mv(yi = lnRR_Sa, V = VCV_S, mod = ~Age_stress_exposure-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat)
summary(mod_S4) ##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -12.4083 24.8166 38.8166 56.1579 40.2166
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0048 0.0694 30 no Study_ID
## sigma^2.2 0.0392 0.1979 92 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 88) = 881.1229, p-val < .0001
##
## Test of Moderators (coefficients 1:4):
## F(df1 = 4, df2 = 88) = 4.3703, p-val = 0.0029
##
## Model Results:
##
## estimate se tval df pval
## Age_stress_exposureAdolescent 0.0074 0.1159 0.0641 88 0.9490
## Age_stress_exposureAdult -0.2279 0.0622 -3.6664 88 0.0004
## Age_stress_exposureEarly postnatal -0.0561 0.0435 -1.2891 88 0.2007
## Age_stress_exposurePrenatal -0.1145 0.0743 -1.5404 88 0.1271
## ci.lb ci.ub
## Age_stress_exposureAdolescent -0.2228 0.2377
## Age_stress_exposureAdult -0.3514 -0.1044 ***
## Age_stress_exposureEarly postnatal -0.1425 0.0304
## Age_stress_exposurePrenatal -0.2621 0.0332
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_S4) ## R2_marginal R2_coditional
## 0.09987307 1.00000000
Age_S <- orchard_plot(mod_S4, mod = "Age_stress_exposure", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 15), # change font sizes
legend.title = element_text(size = 10),
legend.text = element_text(size = 10))
Age_S Type of stress
The type of stressor used
dat5 <- filter(dat, Type_stress_exposure %in% c("Restraint", "Noise", "MS", "Combination"))
VCV_S3 <- impute_covariance_matrix(vi = dat5$lnRRV_S, cluster = dat5$Study_ID, r = 0.5)
mod_S4 <- rma.mv(yi = lnRR_Sa, V = VCV_S3, mod = ~Type_stress_exposure-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat5)
summary(mod_S4) ##
## Multivariate Meta-Analysis Model (k = 85; method: REML)
##
## logLik Deviance AIC BIC AICc
## -11.4138 22.8276 36.8276 53.5887 38.3618
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0115 0.1071 25 no Study_ID
## sigma^2.2 0.0401 0.2002 85 no ES_ID
## sigma^2.3 0.0000 0.0000 4 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 81) = 897.4553, p-val < .0001
##
## Test of Moderators (coefficients 1:4):
## F(df1 = 4, df2 = 81) = 2.8767, p-val = 0.0278
##
## Model Results:
##
## estimate se tval df pval ci.lb
## Type_stress_exposureCombination -0.0500 0.0892 -0.5605 81 0.5767 -0.2274
## Type_stress_exposureMS -0.0539 0.0560 -0.9630 81 0.3384 -0.1653
## Type_stress_exposureNoise -0.1203 0.1036 -1.1608 81 0.2491 -0.3265
## Type_stress_exposureRestraint -0.2101 0.0704 -2.9863 81 0.0037 -0.3501
## ci.ub
## Type_stress_exposureCombination 0.1274
## Type_stress_exposureMS 0.0575
## Type_stress_exposureNoise 0.0859
## Type_stress_exposureRestraint -0.0701 **
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_S4)## R2_marginal R2_coditional
## 0.07029554 1.00000000
Stressor<- orchard_plot(mod_S4, mod = "Type_stress_exposure", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 15), # change font sizes
legend.title = element_text(size = 10),
legend.text = element_text(size = 10))
StressorStess duration
Was the stress acute or chronic?
dat6 <- filter(dat, Stress_duration %in% c("Chronic", "Acute"))
VCV_S4 <- impute_covariance_matrix(vi = dat6$lnRRV_S, cluster = dat6$Study_ID, r = 0.5)
mod_S5 <-rma.mv(yi = lnRR_Sa, V = VCV_S4, mod = ~Stress_duration-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat6)
summary(mod_S5) ##
## Multivariate Meta-Analysis Model (k = 89; method: REML)
##
## logLik Deviance AIC BIC AICc
## -13.7898 27.5796 37.5796 49.9092 38.3204
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0062 0.0786 29 no Study_ID
## sigma^2.2 0.0391 0.1978 89 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 87) = 915.9393, p-val < .0001
##
## Test of Moderators (coefficients 1:2):
## F(df1 = 2, df2 = 87) = 6.6661, p-val = 0.0020
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## Stress_durationAcute 0.0135 0.0659 0.2045 87 0.8384 -0.1176 0.1445
## Stress_durationChronic -0.1360 0.0373 -3.6456 87 0.0005 -0.2101 -0.0618
##
## Stress_durationAcute
## Stress_durationChronic ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_S5) ## R2_marginal R2_coditional
## 0.08245896 1.00000000
Duration_S <- orchard_plot(mod_S5, mod = "Stress_duration", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 15), # change font sizes
legend.title = element_text(size = 10),
legend.text = element_text(size = 10))
Duration_S Multi-moderator model
The akaike weights for the top set of models with AIC < 6
#selecting moderator levels with k >=5
dat_Sfm <- dat %>%
filter(Type_assay %in% c("Recognition", "Habituation", "Conditioning"),
Type_reinforcement %in% c("Appetitive", "Aversive", "Not applicable"),
Type_stress_exposure %in% c("Restraint", "Noise", "MS", "Combination"),
Stress_duration %in% c("Chronic", "Acute"))
VCV_Sfm <- impute_covariance_matrix(vi = dat_Sfm$lnRRV_E, cluster = dat_Sfm$Study_ID, r = 0.5)
mod_Sfm <- rma.mv(yi = lnRR_Sa, V = VCV_Sfm, mod = ~ Type_assay -1 + Learning_vs_memory + Type_reinforcement + Type_stress_exposure + Age_stress_exposure + Stress_duration, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat_Sfm)
#summary(mod_Sfm)
#r2_ml(mod_Sfm)
res_Sfm <- dredge(mod_Sfm, trace=2)
saveRDS(res_Sfm, file = here("Rdata", "res_Sfm.rds"))
# also saving the full model and data
saveRDS(mod_Sfm, file = here("Rdata", "mod_Sfm.rds"))
saveRDS(dat_Sfm, file = here("Rdata", "dat_Sfm.rds"))dat_Sfm <- readRDS(file = here("Rdata", "dat_Sfm.rds"))
mod_Sfm <- readRDS(file = here("Rdata", "mod_Sfm.rds"))
res_Sfm <- readRDS(file = here("Rdata", "res_Sfm.rds"))
res_Sfm2<- subset(res_Sfm, delta <= 6, recalc.weights=FALSE)
importance(res_Sfm2) ## Type_assay Stress_duration Type_stress_exposure
## Sum of weights: 0.91 0.88 0.22
## N containing models: 7 6 2
## Learning_vs_memory Age_stress_exposure Type_reinforcement
## Sum of weights: 0.16 0.04 0.04
## N containing models: 2 1 1
Publication bias & sensitivity analysis
Publication bias
# funnel plot
Funnel_S <- funnel(mod_Sfm, xlab = "lnRR", ylab = "Standard Error")Funnel_S
#calculating inv effective sample size for use in full meta-regression
dat_Sfm$sqrt_inv_e_n <- with(dat_Sfm, sqrt(1/CC_n + 1/EC_n + 1/ES_n + 1/CS_n))
#time lag bias and eggers regression
dat_Sfm$c_Year_published <- as.vector(scale(dat_Sfm$Year_published, scale = F))
PB_MR_S<- rma.mv(lnRR_Sa, lnRRV_S, mods = ~1 + sqrt_inv_e_n + c_Year_published + Type_assay +Learning_vs_memory + Type_reinforcement + Type_stress_exposure + Age_stress_exposure, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
method = "REML",
test = "t",
data = dat_Sfm,
control=list(optimizer="optim", optmethod="Nelder-Mead"))
estimates_PB_MR_S<- estimates.CI(PB_MR_S)
#estimates_PB_MR_SLeave-one-out sensitivity analysis
dat$Study_ID <- as.factor(dat$Study_ID)
LeaveOneOut_effectsize <- list()
for(i in 1:length(levels(dat$Study_ID))){
LeaveOneOut_effectsize[[i]] <- rma.mv(yi = lnRR_Sa, V = lnRRV_S,
random = list(~1 | Study_ID,~1| ES_ID, ~1 | Strain),
method = "REML", data = dat[dat$Study_ID
!= levels(dat$Study_ID)[i], ])}
# writing function for extracting est, ci.lb, and ci.ub from all models
est.func <- function(mod_E0){
df <- data.frame(est = mod_E0$b, lower = mod_E0$ci.lb, upper = mod_E0$ci.ub)
return(df)
}
#using dplyr to form data frame
MA_CVR_S <- lapply(LeaveOneOut_effectsize, function(x) est.func(x))%>% bind_rows %>% mutate(left_out = levels(dat$Study_ID))
saveRDS(MA_CVR_S,file = here("Rdata", "MA_CVR_S.rds"))MA_CVR_S <- readRDS(file = here("Rdata", "MA_CVR_S.rds"))
#telling ggplot to stop reordering factors
MA_CVR_S$left_out<- as.factor(MA_CVR_S$left_out)
MA_CVR_S$left_out<-factor(MA_CVR_S$left_out, levels = MA_CVR_S$left_out)
#plotting
leaveoneout_S <- ggplot(MA_CVR_S) +
geom_hline(yintercept = 0, lty = 2, lwd = 1) +
geom_hline(yintercept = mod_S0$ci.lb, lty = 3, lwd = 0.75, colour = "black") +
geom_hline(yintercept = mod_S0$b, lty = 1, lwd = 0.75, colour = "black") +
geom_hline(yintercept = mod_S0$ci.ub, lty = 3, lwd = 0.75, colour = "black") +
geom_pointrange(aes(x = left_out, y = est, ymin = lower, ymax = upper)) +
xlab("Study left out") +
ylab("lnRR, 95% CI") +
coord_flip() +
theme(panel.grid.minor = element_blank())+
theme_bw() + theme(panel.grid.major = element_blank()) +
theme(panel.grid.minor.x = element_blank() ) +
theme(axis.text.y = element_text(size = 6))
leaveoneout_Sdat$Study_ID <- as.integer(dat$Study_ID)Interaction of stress and EE
Meta-analysis
mod_ES0 <- rma.mv(yi = lnRR_ESa, V = VCV_ES, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat)
summary(mod_ES0) ##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -40.8178 81.6355 89.6355 99.6790 90.1006
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0316 0.1777 30 no Study_ID
## sigma^2.2 0.0229 0.1513 92 no ES_ID
## sigma^2.3 0.0030 0.0544 6 no Strain
##
## Test for Heterogeneity:
## Q(df = 91) = 303.2179, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## 0.1229 0.0596 2.0605 91 0.0422 0.0044 0.2414 *
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
i2_ml(mod_ES0) ## I2_total I2_Study_ID I2_ES_ID I2_Strain
## 0.81703063 0.44913873 0.32576306 0.04212884
orchard_plot(mod_ES0, mod = "Int", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 24), # change font sizes
legend.title = element_text(size = 15),
legend.text = element_text(size = 13)) Meta-regression: uni-moderator
Learning vs Memory
Was the response learning or memory?
mod_ES1 <- rma.mv(yi = lnRR_ESa, V = VCV_ES, mod = ~Learning_vs_memory-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat)
summary(mod_ES1) ##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -40.4769 80.9539 90.9539 103.4529 91.6682
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0292 0.1708 30 no Study_ID
## sigma^2.2 0.0232 0.1524 92 no ES_ID
## sigma^2.3 0.0034 0.0582 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 90) = 299.1854, p-val < .0001
##
## Test of Moderators (coefficients 1:2):
## F(df1 = 2, df2 = 90) = 2.9219, p-val = 0.0590
##
## Model Results:
##
## estimate se tval df pval ci.lb
## Learning_vs_memoryLearning 0.1744 0.0722 2.4166 90 0.0177 0.0310
## Learning_vs_memoryMemory 0.1057 0.0619 1.7065 90 0.0914 -0.0174
## ci.ub
## Learning_vs_memoryLearning 0.3178 *
## Learning_vs_memoryMemory 0.2288 .
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_ES1) ## R2_marginal R2_coditional
## 0.0197648 0.9405080
LvsM_ES <- orchard_plot(mod_ES1, mod = "Learning_vs_memory", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
geom_point(aes(fill = name), size = 3, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
axis.text.x = element_text(size = 10), # change font sizes
axis.text.y = element_text(size = 10),
legend.title = element_text(size = 7),
legend.text = element_text(size = 7))
LvsM_ES Type of assay
The broad category of the type of assay used to measure learning or memory
VCV_ES1 <- impute_covariance_matrix(vi = dat1$lnRRV_ES, cluster = dat$Study_ID, r = 0.5)
mod_ES2 <- rma.mv(yi = lnRR_ESa, V = VCV_ES1, mod = ~Type_assay-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat1)
summary(mod_ES2)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -39.0874 78.1748 90.1748 105.1066 91.1992
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0370 0.1924 30 no Study_ID
## sigma^2.2 0.0192 0.1386 92 no ES_ID
## sigma^2.3 0.0018 0.0422 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 89) = 293.9385, p-val < .0001
##
## Test of Moderators (coefficients 1:3):
## F(df1 = 3, df2 = 89) = 3.1062, p-val = 0.0305
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## Type_assayConditioning 0.1525 0.0589 2.5877 89 0.0113 0.0354 0.2696
## Type_assayHabituation 0.1990 0.1415 1.4070 89 0.1629 -0.0820 0.4801
## Type_assayRecognition -0.0048 0.0800 -0.0606 89 0.9518 -0.1637 0.1541
##
## Type_assayConditioning *
## Type_assayHabituation
## Type_assayRecognition
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_ES2) ## R2_marginal R2_coditional
## 0.05775809 0.97105550
Learning_ES <- orchard_plot(mod_ES2, mod = "Type_assay", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
geom_point(aes(fill = name), size = 3, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
axis.text.x = element_text(size = 10), # change font sizes
axis.text.y = element_text(size = 10),
legend.title = element_text(size = 7),
legend.text = element_text(size = 7))
Learning_ESType of reinforcement
If conditioning was used, was aversive or appetitive reinforcement used?
VCV_ES2 <- impute_covariance_matrix(vi = dat2$lnRRV_ES, cluster = dat2$Study_ID, r = 0.5)
mod_ES3 <- rma.mv(yi = lnRR_ESa, V = VCV_ES2, mod = ~ Type_reinforcement-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat2)
summary(mod_ES3)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -39.0604 78.1208 90.1208 105.0526 91.1452
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0382 0.1954 30 no Study_ID
## sigma^2.2 0.0189 0.1377 92 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 89) = 293.4724, p-val < .0001
##
## Test of Moderators (coefficients 1:3):
## F(df1 = 3, df2 = 89) = 3.2547, p-val = 0.0254
##
## Model Results:
##
## estimate se tval df pval ci.lb
## Type_reinforcementAppetitive 0.1007 0.1075 0.9366 89 0.3515 -0.1129
## Type_reinforcementAversive 0.1573 0.0569 2.7618 89 0.0070 0.0441
## Type_reinforcementNot applicable 0.0147 0.0702 0.2101 89 0.8341 -0.1247
## ci.ub
## Type_reinforcementAppetitive 0.3143
## Type_reinforcementAversive 0.2704 **
## Type_reinforcementNot applicable 0.1541
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_ES3) ## R2_marginal R2_coditional
## 0.0586952 1.0000000
Reinforcement_ES <- orchard_plot(mod_ES3, mod = "Type_reinforcement", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
geom_point(aes(fill = name), size = 3, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
axis.text.x = element_text(size = 10), # change font sizes
axis.text.y = element_text(size = 10),
legend.title = element_text(size = 7),
legend.text = element_text(size = 7))
Reinforcement_ES Age of enrichment
The age when individuals were exposed to environmental enrichment
VCV_ES3 <- impute_covariance_matrix(vi = dat3$lnRRV_ES, cluster = dat3$Study_ID, r = 0.5)
mod_ES4 <- rma.mv(yi = lnRR_ESa, V = VCV_ES3, mod = ~Age_EE_exposure-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat3)
summary(mod_ES4) ##
## Multivariate Meta-Analysis Model (k = 88; method: REML)
##
## logLik Deviance AIC BIC AICc
## -36.6407 73.2813 83.2813 95.5531 84.0313
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0332 0.1822 29 no Study_ID
## sigma^2.2 0.0207 0.1439 88 no ES_ID
## sigma^2.3 0.0007 0.0265 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 86) = 288.6493, p-val < .0001
##
## Test of Moderators (coefficients 1:2):
## F(df1 = 2, df2 = 86) = 3.1308, p-val = 0.0487
##
## Model Results:
##
## estimate se tval df pval ci.lb
## Age_EE_exposureAdolescent 0.1284 0.0583 2.2006 86 0.0304 0.0124
## Age_EE_exposureAdult 0.1247 0.0921 1.3538 86 0.1793 -0.0584
## ci.ub
## Age_EE_exposureAdolescent 0.2444 *
## Age_EE_exposureAdult 0.3077
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_ES4) ## R2_marginal R2_coditional
## 0.0000400928 0.9871095822
Age_enrichment_ES <- orchard_plot(mod_ES4, mod = "Age_EE_exposure", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
geom_point(aes(fill = name), size = 3, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
axis.text.x = element_text(size = 10), # change font sizes
axis.text.y = element_text(size = 10),
legend.title = element_text(size = 7),
legend.text = element_text(size = 7))
Age_enrichment_ESExercise enrichment
Did enrichment involve the addition of apparatus for voluntary exercise (e.g., a wheel or treadmill)?
mod_ES5<- rma.mv(yi = lnRR_ESa, V = VCV_ES, mod = ~EE_exercise-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat)
summary(mod_ES5)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -40.4968 80.9935 90.9935 103.4926 91.7078
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0325 0.1803 30 no Study_ID
## sigma^2.2 0.0230 0.1517 92 no ES_ID
## sigma^2.3 0.0065 0.0805 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 90) = 297.3270, p-val < .0001
##
## Test of Moderators (coefficients 1:2):
## F(df1 = 2, df2 = 90) = 1.8401, p-val = 0.1647
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## EE_exerciseExercise 0.1051 0.0780 1.3474 90 0.1812 -0.0499 0.2602
## EE_exerciseNo exercise 0.1687 0.0967 1.7448 90 0.0844 -0.0234 0.3609
##
## EE_exerciseExercise
## EE_exerciseNo exercise .
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_ES5) ## R2_marginal R2_coditional
## 0.01474459 0.89712469
Exercise_ES <- orchard_plot(mod_ES5, mod = "EE_exercise", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
geom_point(aes(fill = name), size = 3, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
axis.text.x = element_text(size = 10), # change font sizes
axis.text.y = element_text(size = 10),
legend.title = element_text(size = 7),
legend.text = element_text(size = 7))
Exercise_ESAge of stress
The age when individuals were exposed to stress
mod_ES7 <-rma.mv(yi = lnRR_ESa, V = VCV_ES, mod = ~Age_stress_exposure-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat)
summary(mod_ES7) ##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -39.0917 78.1834 92.1834 109.5247 93.5834
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0295 0.1718 30 no Study_ID
## sigma^2.2 0.0232 0.1522 92 no ES_ID
## sigma^2.3 0.0091 0.0954 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 88) = 286.4225, p-val < .0001
##
## Test of Moderators (coefficients 1:4):
## F(df1 = 4, df2 = 88) = 1.5932, p-val = 0.1832
##
## Model Results:
##
## estimate se tval df pval
## Age_stress_exposureAdolescent -0.0137 0.1762 -0.0775 88 0.9384
## Age_stress_exposureAdult 0.1677 0.1140 1.4708 88 0.1449
## Age_stress_exposureEarly postnatal 0.1067 0.0920 1.1602 88 0.2491
## Age_stress_exposurePrenatal 0.3179 0.1311 2.4254 88 0.0173
## ci.lb ci.ub
## Age_stress_exposureAdolescent -0.3639 0.3366
## Age_stress_exposureAdult -0.0589 0.3942
## Age_stress_exposureEarly postnatal -0.0761 0.2896
## Age_stress_exposurePrenatal 0.0574 0.5784 *
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_ES7) ## R2_marginal R2_coditional
## 0.09276574 0.86630830
Age_stress_ES<-orchard_plot(mod_ES7, mod = "Age_stress_exposure", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
geom_point(aes(fill = name), size = 3, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
axis.text.x = element_text(size = 10), # change font sizes
axis.text.y = element_text(size = 10),
legend.title = element_text(size = 7),
legend.text = element_text(size = 7))
Age_stress_ESType of stress
The type of stressor used
VCV_ES5 <- impute_covariance_matrix(vi = dat5$lnRRV_ES, cluster = dat5$Study_ID, r = 0.5)
mod_ES8 <- rma.mv(yi = lnRR_ESa, V = VCV_ES5, mod = ~Type_stress_exposure-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat5)
summary(mod_ES8)##
## Multivariate Meta-Analysis Model (k = 85; method: REML)
##
## logLik Deviance AIC BIC AICc
## -34.4046 68.8091 82.8091 99.5703 84.3434
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0426 0.2064 25 no Study_ID
## sigma^2.2 0.0232 0.1524 85 no ES_ID
## sigma^2.3 0.0104 0.1021 4 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 81) = 281.9708, p-val < .0001
##
## Test of Moderators (coefficients 1:4):
## F(df1 = 4, df2 = 81) = 0.5137, p-val = 0.7258
##
## Model Results:
##
## estimate se tval df pval ci.lb
## Type_stress_exposureCombination 0.1111 0.1458 0.7618 81 0.4484 -0.1790
## Type_stress_exposureMS 0.1185 0.1099 1.0784 81 0.2840 -0.1001
## Type_stress_exposureNoise 0.1651 0.1795 0.9198 81 0.3604 -0.1920
## Type_stress_exposureRestraint 0.1374 0.1252 1.0978 81 0.2755 -0.1116
## ci.ub
## Type_stress_exposureCombination 0.4011
## Type_stress_exposureMS 0.3370
## Type_stress_exposureNoise 0.5221
## Type_stress_exposureRestraint 0.3865
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_ES8)## R2_marginal R2_coditional
## 0.004455703 0.863910284
Stressor_ES <- orchard_plot(mod_ES8, mod = "Type_stress_exposure", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
geom_point(aes(fill = name), size = 3, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
axis.text.x = element_text(size = 10), # change font sizes
axis.text.y = element_text(size = 10),
legend.title = element_text(size = 7),
legend.text = element_text(size = 7))
Stressor_ES Stress duration
Was the stress acute or chronic?
VCV_ES6 <- impute_covariance_matrix(vi = dat6$lnRRV_ES, cluster = dat6$Study_ID, r = 0.5)
mod_ES9 <-rma.mv(yi = lnRR_ESa, V = VCV_ES6, mod = ~Stress_duration-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat6)
summary(mod_ES9) ##
## Multivariate Meta-Analysis Model (k = 89; method: REML)
##
## logLik Deviance AIC BIC AICc
## -35.9010 71.8020 81.8020 94.1315 82.5427
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0133 0.1155 29 no Study_ID
## sigma^2.2 0.0260 0.1611 89 no ES_ID
## sigma^2.3 0.0080 0.0895 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 87) = 278.4877, p-val < .0001
##
## Test of Moderators (coefficients 1:2):
## F(df1 = 2, df2 = 87) = 4.3877, p-val = 0.0153
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## Stress_durationAcute -0.0367 0.0930 -0.3946 87 0.6941 -0.2216 0.1482
## Stress_durationChronic 0.1854 0.0732 2.5347 87 0.0130 0.0400 0.3308
##
## Stress_durationAcute
## Stress_durationChronic *
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_ES9) ## R2_marginal R2_coditional
## 0.1598311 0.8575918
Duration_ES<- orchard_plot(mod_ES9, mod = "Stress_duration", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
geom_point(aes(fill = name), size = 3, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
axis.text.x = element_text(size = 10), # change font sizes
axis.text.y = element_text(size = 10),
legend.title = element_text(size = 7),
legend.text = element_text(size = 7))
Duration_ESOrder to treatment exposure
The order in which individuals were exposed to enrichment and stress
mod_ES10 <- rma.mv(yi = lnRR_ESa, V = VCV_ES, mod = ~Exposure_order -1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat)
summary(mod_ES10)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -39.0408 78.0817 90.0817 105.0135 91.1061
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0316 0.1777 30 no Study_ID
## sigma^2.2 0.0227 0.1507 92 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 89) = 292.2561, p-val < .0001
##
## Test of Moderators (coefficients 1:3):
## F(df1 = 3, df2 = 89) = 3.1546, p-val = 0.0287
##
## Model Results:
##
## estimate se tval df pval ci.lb
## Exposure_orderConcurrently 0.1492 0.1208 1.2351 89 0.2201 -0.0909
## Exposure_orderEnrichment first -0.1782 0.1659 -1.0744 89 0.2856 -0.5079
## Exposure_orderStress first 0.1370 0.0526 2.6046 89 0.0108 0.0325
## ci.ub
## Exposure_orderConcurrently 0.3893
## Exposure_orderEnrichment first 0.1514
## Exposure_orderStress first 0.2414 *
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_ES10)## R2_marginal R2_coditional
## 0.1027207 1.0000000
Order_ES <- orchard_plot(mod_ES10, mod = "Exposure_order", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
geom_point(aes(fill = name), size = 3, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
axis.text.x = element_text(size = 10), # change font sizes
axis.text.y = element_text(size = 10),
legend.title = element_text(size = 7),
legend.text = element_text(size = 7))
Order_ES Multi-moderator model
The akaike weights for the top set of models with AIC < 6
dat_ESfm <- dat %>%
filter(Type_assay %in% c("Recognition", "Habituation", "Conditioning"),
Type_reinforcement %in% c("Appetitive", "Aversive", "Not applicable"),
EE_social %in% c("Social", "Non-social"),
Age_EE_exposure %in% c("Adult", "Adolescent"),
Type_stress_exposure %in% c("Restraint", "Noise", "MS", "Combination"),
Stress_duration %in% c("Chronic", "Acute"),
Age_stress_exposure %in% c("Prenatal", "Early postnatal", "Adult"))
VCV_ESfm <- impute_covariance_matrix(vi = dat_ESfm$lnRRV_ES, cluster = dat_ESfm$Study_ID, r = 0.5)
mod_ESfm <- rma.mv(yi = lnRR_Sa, V = VCV_ESfm, mod = ~Type_assay-1 + Learning_vs_memory + Type_reinforcement + EE_social + EE_exercise + Age_EE_exposure + Type_stress_exposure + Age_stress_exposure + Stress_duration + Exposure_order, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat_ESfm)
#summary(mod_ESfm)
#r2_ml(mod_ESfm)
res_ESfm <- dredge(mod_ESfm, trace=2)
saveRDS(res_ESfm, file = here("Rdata", "res_ESfm.rds"))
# also saving the full model and data
saveRDS(mod_ESfm, file = here("Rdata", "mod_ESfm.rds"))
saveRDS(dat_ESfm, file = here("Rdata", "dat_ESfm.rds"))dat_ESfm <- readRDS(file = here("Rdata", "dat_ESfm.rds"))
mod_ESfm <- readRDS(file = here("Rdata", "mod_ESfm.rds"))
res_ESfm <- readRDS(file = here("Rdata", "res_ESfm.rds"))
res_ESfm2<- subset(res_ESfm, delta <= 6, recalc.weights=FALSE)
importance(res_ESfm2) ## Type_assay Stress_duration Age_EE_exposure EE_exercise
## Sum of weights: 0.65 0.63 0.26 0.11
## N containing models: 19 18 7 6
## Learning_vs_memory EE_social Age_stress_exposure
## Sum of weights: 0.08 0.08 0.05
## N containing models: 4 4 2
## Type_stress_exposure Exposure_order
## Sum of weights: 0.03 0.02
## N containing models: 2 1
Publication bias & sensitivity analysis
Publication bias
Funnel_ES<-funnel(mod_ESfm, xlab = "lnRR", ylab = "Standard Error")Funnel_ES
#year published was scaled previously under stress PB
dat_ESfm$sqrt_inv_e_n <- with(dat_ESfm, sqrt(1/CC_n + 1/EC_n + 1/ES_n + 1/CS_n))
PB_MR_ES<- rma.mv(lnRR_ESa, lnRRV_ES, mods = ~1 + sqrt_inv_e_n + Year_published + Learning_vs_memory + Type_assay + Type_reinforcement + EE_social + EE_exercise + Age_stress_exposure, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain), method = "REML", test = "t",
data = dat_ESfm)
estimates_PB_MR_ES<- estimates.CI(PB_MR_ES)
#estimates_PB_MR_ESLeave-one-out analysis
dat$Study_ID <- as.factor(dat$Study_ID)
LeaveOneOut_effectsize <- list()
for(i in 1:length(levels(dat$Study_ID))){
LeaveOneOut_effectsize[[i]] <- rma.mv(yi = lnRR_Ea, V = lnRRV_E,
random = list(~1 | Study_ID,~1| ES_ID, ~1 | Strain),
method = "REML", data = dat[dat$Study_ID
!= levels(dat$Study_ID)[i], ])}
# writing function for extracting est, ci.lb, and ci.ub from all models
est.func <- function(mod_E0){
df <- data.frame(est = mod_E0$b, lower = mod_E0$ci.lb, upper = mod_E0$ci.ub)
return(df)
}
#using dplyr to form data frame
MA_CVR_ES <- lapply(LeaveOneOut_effectsize, function(x) est.func(x))%>% bind_rows %>% mutate(left_out = levels(dat$Study_ID))
saveRDS(MA_CVR_ES, ,file = here("Rdata", "MA_CVR_ES.rds"))MA_CVR_ES<- readRDS(here("Rdata", "MA_CVR_ES.rds"))
#telling ggplot to stop reordering factors
MA_CVR_ES$left_out<- as.factor(MA_CVR_ES$left_out)
MA_CVR_ES$left_out<-factor(MA_CVR_ES$left_out, levels = MA_CVR_ES$left_out)
#plotting
leaveoneout_ES <- ggplot(MA_CVR_ES) +
geom_hline(yintercept = 0, lty = 2, lwd = 1) +
geom_hline(yintercept = mod_E0$ci.lb, lty = 3, lwd = 0.75, colour = "black") +
geom_hline(yintercept = mod_E0$b, lty = 1, lwd = 0.75, colour = "black") +
geom_hline(yintercept = mod_E0$ci.ub, lty = 3, lwd = 0.75, colour = "black") +
geom_pointrange(aes(x = left_out, y = est, ymin = lower, ymax = upper)) +
xlab("Study left out") +
ylab("lnRR, 95% CI") +
coord_flip() +
theme(panel.grid.minor = element_blank())+
theme_bw() + theme(panel.grid.major = element_blank()) +
theme(panel.grid.minor.x = element_blank() ) +
theme(axis.text.y = element_text(size = 6))
leaveoneout_ESdat$Study_ID <- as.integer(dat$Study_ID)Combined orchard plot
mod_list1 <- list(mod_E0, mod_S0, mod_ES0)
mod_res1 <- lapply(mod_list1, function(x) mod_results(x, mod = "Int"))
merged1 <- submerge(mod_res1[[3]], mod_res1[[2]], mod_res1[[1]], mix = T)
merged1$mod_table$name <- factor(merged1$mod_table$name, levels = c("Intrcpt1",
"Intrcpt2", "Intrcpt3"),
labels = rev(c("Enrichment ME", "Stress ME", "Interaction")))
merged1$data$moderator <- factor(merged1$data$moderator, levels = c("Intrcpt1",
"Intrcpt2", "Intrcpt3"),
labels = rev(c("Enrichment ME", "Stress ME", "Interaction")))
orchard1<- orchard_plot(merged1, mod = "Int", xlab = "lnRR", angle = 0) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
xlim(-2,4.5) +
geom_point(aes(fill = name), size = 4, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling +
scale_colour_manual(values = c("#00AEEF","#00A651","#ED1C24"))+ # change colours
scale_fill_manual(values=c("#00AEEF","#00A651","#ED1C24"))+
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 15), # change font sizes
legend.title = element_text(size = 10),
legend.text = element_text(size = 10))
orchard1‘Pairwise’ effect sizes
Enrichment relative to control
VCV_E20 <- impute_covariance_matrix(vi = dat$lnRRV_E2, cluster = dat$Study_ID, r = 0.5)
mod_E20 <- rma.mv(yi = lnRR_E2a, V = VCV_E20, random = list(~1|Study_ID,
~ 1|Strain,
~1|ES_ID),
test = "t",
data = dat,
control=list(optimizer="optim", optmethod="Nelder-Mead"))
summary(mod_E20) ##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -7.3235 14.6470 22.6470 32.6904 23.1121
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0037 0.0611 30 no Study_ID
## sigma^2.2 0.0000 0.0000 6 no Strain
## sigma^2.3 0.0281 0.1675 92 no ES_ID
##
## Test for Heterogeneity:
## Q(df = 91) = 475.8327, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## 0.1066 0.0291 3.6655 91 0.0004 0.0489 0.1644 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
i2_ml(mod_E20) ## I2_total I2_Study_ID I2_Strain I2_ES_ID
## 8.610789e-01 1.011724e-01 2.607945e-09 7.599065e-01
orchard_plot(mod_E20, mod = "Int", xlab = "lnRR")Stress relative to control
VCV_S20 <- impute_covariance_matrix(vi = dat$lnRRV_S2, cluster = dat$Study_ID, r = 0.5)
mod_S20 <- rma.mv(yi = lnRR_S2a, V = VCV_S20, random = list(~1|Study_ID,
~ 1|Strain,
~1|ES_ID),
test = "t",
data = dat)
summary(mod_S20) ##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -52.3561 104.7122 112.7122 122.7557 113.1773
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0323 0.1797 30 no Study_ID
## sigma^2.2 0.0000 0.0000 6 no Strain
## sigma^2.3 0.0798 0.2824 92 no ES_ID
##
## Test for Heterogeneity:
## Q(df = 91) = 1003.0694, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## -0.1846 0.0522 -3.5360 91 0.0006 -0.2882 -0.0809 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
i2_ml(mod_S20) ## I2_total I2_Study_ID I2_Strain I2_ES_ID
## 9.489604e-01 2.734220e-01 9.879730e-10 6.755384e-01
orchard_plot(mod_S20, mod = "Int", xlab = "lnRR")Enrichment + stress relative to control
VCV_ES20 <- impute_covariance_matrix(vi = dat$lnRRV_ES2, cluster = dat$Study_ID, r = 0.5)
mod_ES20 <- rma.mv(yi = lnRR_ES2a, V = VCV_ES20, random = list(~1|Study_ID,
~ 1|Strain,
~1|ES_ID),
test = "t",
data = dat)
summary(mod_ES20) ##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -10.3625 20.7250 28.7250 38.7684 29.1901
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0039 0.0625 30 no Study_ID
## sigma^2.2 0.0014 0.0377 6 no Strain
## sigma^2.3 0.0227 0.1508 92 no ES_ID
##
## Test for Heterogeneity:
## Q(df = 91) = 402.2656, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## 0.0801 0.0389 2.0594 91 0.0423 0.0028 0.1573 *
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
i2_ml(mod_ES20) ## I2_total I2_Study_ID I2_Strain I2_ES_ID
## 0.81513970 0.11355557 0.04132363 0.66026050
orchard_plot(mod_ES20, mod = "Int", xlab = "lnRR")Enrichment + stress relative to stress
VCV_E30 <- impute_covariance_matrix(vi = dat$lnRRV_E3, cluster = dat$Study_ID, r = 0.5)
mod_E30 <- rma.mv(yi = lnRR_E3a, V = VCV_E30, random = list(~1|Study_ID,
~ 1|Strain,
~1|ES_ID),
test = "t",
data = dat)
summary(mod_E30) ##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -46.3447 92.6895 100.6895 110.7329 101.1546
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0235 0.1532 30 no Study_ID
## sigma^2.2 0.0263 0.1623 6 no Strain
## sigma^2.3 0.0543 0.2331 92 no ES_ID
##
## Test for Heterogeneity:
## Q(df = 91) = 790.0249, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## 0.2465 0.1011 2.4389 91 0.0167 0.0457 0.4472 *
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
i2_ml(mod_E30) ## I2_total I2_Study_ID I2_Strain I2_ES_ID
## 0.9456920 0.2132063 0.2391809 0.4933048
orchard_plot(mod_E30, mod = "Int", xlab = "lnRR")Enrichment + stress relative to enrichment
VCV_S30 <- impute_covariance_matrix(vi = dat$lnRRV_S3, cluster = dat$Study_ID, r = 0.5)
mod_S30 <- rma.mv(yi = lnRR_S3a, V = VCV_S30, random = list(~1|Study_ID,
~ 1|Strain,
~1|ES_ID),
test = "t",
data = dat,
control=list(optimizer="optim", optmethod="Nelder-Mead"))
summary(mod_S30) ##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -6.8788 13.7576 21.7576 31.8011 22.2228
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0000 0.0000 30 no Study_ID
## sigma^2.2 0.0009 0.0292 6 no Strain
## sigma^2.3 0.0276 0.1662 92 no ES_ID
##
## Test for Heterogeneity:
## Q(df = 91) = 540.3522, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## -0.0087 0.0342 -0.2552 91 0.7992 -0.0766 0.0592
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
i2_ml(mod_S30) ## I2_total I2_Study_ID I2_Strain I2_ES_ID
## 8.415428e-01 5.192043e-09 2.515933e-02 8.163835e-01
orchard_plot(mod_S30, mod = "Int", xlab = "lnRR")mod_list2 <- list(mod_S30, mod_E30, mod_ES20, mod_S20, mod_E20) #rearranged the order so that it matches intext results
mod_res2 <- lapply(mod_list2, function(x) mod_results(x, mod = "Int"))
merged2 <- submerge(mod_res2[[1]], mod_res2[[2]], mod_res2[[3]], mod_res2[[4]], mod_res2[[5]], mix = T)
merged2$mod_table$name <- factor(merged2$mod_table$name, levels = c("Intrcpt1",
"Intrcpt2", "Intrcpt3", "Intrcpt4", "Intrcpt5"),
labels = rev(c("EC/CC", "CS/CC", "ES/CC", "ES/CS", "ES/EC")))
merged2$data$moderator <- factor(merged2$data$moderator, levels = c("Intrcpt1",
"Intrcpt2", "Intrcpt3", "Intrcpt4", "Intrcpt5"),
labels = rev(c("EC/CC", "CS/CC", "ES/CC", "ES/CS", "ES/EC")))
orchard2 <- orchard_plot(merged2, mod = "Int", xlab = "lnRR", angle = 0) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
xlim(-2,4.5) +
geom_point(aes(fill = name), size = 4, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
scale_colour_manual(values = c("#7B81BE","#D7DF23","#F37158","#75CBF2","#97D2B4"))+ # change colours
scale_fill_manual(values=c("#7B81BE","#D7DF23","#F37158","#75CBF2","#97D2B4"))+
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 15), # change font sizes
legend.title = element_text(size = 10),
legend.text = element_text(size = 10)) Figures
Panel of ‘focal’ ES and ‘pairwise’ ES orchard plots
p1 <- orchard1 + orchard2 + plot_annotation(tag_levels = 'A')
p1#saved as PDF: 6 x 15 inchesPanel of meta-regressions
Environmental enrichment
#Enrichment
E_mod <- (LvsM_E + Learning_E + Reinforcement_E)/ (Age_E + Exercise_E + Social_E) + plot_annotation(tag_levels = 'A')
E_mod#saved as pdf 10 x 15 inchesStress
S_mod <- (LvsM_S + Learning_S + Reinforcement_S) / (Age_S + Stressor + Duration_S) + plot_annotation(tag_levels = 'A')
S_mod#saved as pdf 10 x 15 inchesInteraction
ES_mod <- plot_grid(LvsM_ES, Learning_ES, Reinforcement_ES, Age_enrichment_ES, Age_stress_ES, Order_ES, Exercise_ES, Social_ES, Stressor_ES, Duration_ES,
labels = "AUTO", ncol = 5)
ES_mod#saved as 10 x 20 inchesPanel of funnel plots
# EE
pdf(NULL)
dev.control(displaylist="enable")
par(mar=c(4,4,0.1,0))
A <- funnel(mod_Sfm, xlab = "Residuals (lnRR)", ylab = "Standard Error",
xlim = c(-2,2),
ylim = c(0,1.05))
A <- recordPlot()
invisible(dev.off())
# Stress
pdf(NULL)
dev.control(displaylist="enable")
par(mar=c(4,4,0.1,0))
B <- funnel(mod_Sfm, xlab = "Residuals (lnRR)", ylab = "Standard Error",
xlim = c(-2,2),
ylim = c(0,1.05))
B <- recordPlot()
invisible(dev.off())
# Interaction
pdf(NULL)
dev.control(displaylist="enable")
par(mar=c(4,4,0.1,0))
C <- funnel(mod_ESfm, xlab = "Residuals (lnRR)", ylab = "Standard Error",
xlim = c(-2,2),
ylim = c(0,1.05))
C <- recordPlot()
invisible(dev.off())
# putting together
ggdraw(A) + ggdraw(B) + ggdraw(C) + plot_annotation(tag_levels = 'A')
#png(file = here("figs", "Fig7_Funnels.png"))
#dev.off()knitr::include_graphics(here("figs", "funnels.png"))Modelling with SMD
Models are of the three ‘focal’ models
Environmental Enrichment
mod_E0a <- rma.mv(yi = SMD_Ea, V = VCV_Ea, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat)
summary(mod_E0a)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -111.6284 223.2568 231.2568 241.3003 231.7220
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.2512 0.5012 30 no Study_ID
## sigma^2.2 0.4247 0.6517 92 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Heterogeneity:
## Q(df = 91) = 673.4722, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## 0.7241 0.1295 5.5895 91 <.0001 0.4668 0.9814 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
i2_ml(mod_E0a)## I2_total I2_Study_ID I2_ES_ID I2_Strain
## 8.673320e-01 3.223202e-01 5.450118e-01 1.749563e-09
Stress
mod_S0a <- rma.mv(yi = SMD_Sa, V = VCV_Sa, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat)
summary(mod_S0a) ##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -120.2817 240.5634 248.5634 258.6068 249.0285
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.3874 0.6224 30 no Study_ID
## sigma^2.2 0.4969 0.7049 92 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Heterogeneity:
## Q(df = 91) = 818.5592, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## -0.4256 0.1498 -2.8403 91 0.0056 -0.7232 -0.1279 **
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
i2_ml(mod_S0a) ## I2_total I2_Study_ID I2_ES_ID I2_Strain
## 8.959540e-01 3.924790e-01 5.034750e-01 1.438454e-09
Interaction
mod_ES0a <- rma.mv(yi = SMD_ESa, V = VCV_ESa, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat)
summary(mod_ES0a)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -126.0571 252.1141 260.1141 270.1576 260.5793
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.4653 0.6821 30 no Study_ID
## sigma^2.2 0.3698 0.6081 92 no ES_ID
## sigma^2.3 0.0000 0.0002 6 no Strain
##
## Test for Heterogeneity:
## Q(df = 91) = 257.4673, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## 0.6880 0.1763 3.9017 91 0.0002 0.3377 1.0382 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
i2_ml(mod_ES0a) ## I2_total I2_Study_ID I2_ES_ID I2_Strain
## 6.657130e-01 3.709364e-01 2.947766e-01 3.351907e-08
Risk of Bias
Percent of studies that used randomisation and used blinding
# randomisation
dat %>% group_by(ROB_randomisation) %>%
summarise(n = n_distinct(Study_ID))## # A tibble: 2 x 2
## ROB_randomisation n
## <chr> <int>
## 1 Unclear 15
## 2 Yes 15
15/30## [1] 0.5
#blinding
dat %>% group_by(ROB_blinding) %>%
summarise(n = n_distinct(Study_ID))## # A tibble: 2 x 2
## ROB_blinding n
## <chr> <int>
## 1 Unclear 24
## 2 Yes 6
6/30## [1] 0.2
Software and package versions
sessionInfo() %>% pander()R version 4.1.0 (2021-05-18)
Platform: x86_64-apple-darwin17.0 (64-bit)
locale: en_AU.UTF-8||en_AU.UTF-8||en_AU.UTF-8||C||en_AU.UTF-8||en_AU.UTF-8
attached base packages: grid, stats, graphics, grDevices, utils, datasets, methods and base
other attached packages: pander(v.0.6.4), gridGraphics(v.0.5-1), png(v.0.1-7), cowplot(v.1.1.1), ggthemr(v.1.1.0), ggalluvial(v.0.12.3), visdat(v.0.5.3), networkD3(v.0.4), GoodmanKruskal(v.0.0.3), patchwork(v.1.1.1), MuMIn(v.1.43.17), orchaRd(v.0.0.0.9000), clubSandwich(v.0.5.3), metafor(v.3.0-2), Matrix(v.1.3-3), here(v.1.0.1), forcats(v.0.5.1), stringr(v.1.4.0), dplyr(v.1.0.6), purrr(v.0.3.4), readr(v.1.4.0), tidyr(v.1.1.3), tibble(v.3.1.2), ggplot2(v.3.3.3) and tidyverse(v.1.3.1)
loaded via a namespace (and not attached): nlme(v.3.1-152), fs(v.1.5.0), lubridate(v.1.7.10), httr(v.1.4.2), rprojroot(v.2.0.2), tools(v.4.1.0), backports(v.1.2.1), utf8(v.1.2.1), R6(v.2.5.0), DBI(v.1.1.1), colorspace(v.2.0-1), withr(v.2.4.2), tidyselect(v.1.1.1), compiler(v.4.1.0), cli(v.2.5.0), rvest(v.1.0.0), pacman(v.0.5.1), xml2(v.1.3.2), sandwich(v.3.0-1), bookdown(v.0.22), scales(v.1.1.1), digest(v.0.6.27), rmarkdown(v.2.8), pkgconfig(v.2.0.3), htmltools(v.0.5.1.1), highr(v.0.9), dbplyr(v.2.1.1), htmlwidgets(v.1.5.3), rlang(v.0.4.11), readxl(v.1.3.1), rstudioapi(v.0.13), generics(v.0.1.0), zoo(v.1.8-9), jsonlite(v.1.7.2), magrittr(v.2.0.1), Rcpp(v.1.0.6), munsell(v.0.5.0), fansi(v.0.4.2), lifecycle(v.1.0.0), stringi(v.1.6.2), yaml(v.2.2.1), mathjaxr(v.1.4-0), crayon(v.1.4.1), lattice(v.0.20-44), haven(v.2.4.1), hms(v.1.1.0), knitr(v.1.33), pillar(v.1.6.1), igraph(v.1.2.6), codetools(v.0.2-18), stats4(v.4.1.0), reprex(v.2.0.0), glue(v.1.4.2), evaluate(v.0.14), modelr(v.0.1.8), vctrs(v.0.3.8), rmdformats(v.1.0.2), cellranger(v.1.1.0), gtable(v.0.3.0), assertthat(v.0.2.1), xfun(v.0.23), broom(v.0.7.6) and ellipsis(v.0.3.2)
Social enrichment
Did enrichment involve more conspecifics (note that we did not include studies that only provided social enrichment but no other form of abiotic enrichment)?